﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "rank")]
    public static unsafe int rank(IntPtr a_ptr, int m, int n)
    {
        double* a = (double*)a_ptr.ToPointer();

        return rank(a, m, n);
    }

    /// <summary>
    /// 矩阵求秩
    /// </summary>
    /// <param name="a">a[m][n]存放m×n阶矩阵A的元素。返回时将被破坏。</param>
    /// <param name="m"></param>
    /// <param name="n"></param>
    /// <returns>返回A的秩。</returns>
    public static unsafe int rank(double* a, int m, int n)
    {
        int i, j, k, nn, si = 0, js = 0, l, ll, u, v;
        double q, d;

        nn = m;
        if (m >= n) nn = n;
        k = 0;
        for (l = 0; l <= nn - 1; l++)
        {
            q = 0.0;
            for (i = l; i <= m - 1; i++)
            {
                for (j = l; j <= n - 1; j++)
                {
                    ll = i * n + j;
                    d = Math.Abs(a[ll]);
                    if (d > q)
                    {
                        q = d;
                        si = i;
                        js = j;
                    }
                }
            }
            //if (Math.Abs(q) <= float.Epsilon)
            if (EQ0(q))
            {
                return (k);
            }
            k = k + 1;
            if (si != l)
            {
                for (j = l; j <= n - 1; j++)
                {
                    u = l * n + j;
                    v = si * n + j;
                    d = a[u]; a[u] = a[v]; a[v] = d;
                }
            }
            if (js != l)
            {
                for (i = l; i <= m - 1; i++)
                {
                    u = i * n + js;
                    v = i * n + l;
                    d = a[u]; a[u] = a[v]; a[v] = d;
                }
            }
            ll = l * n + l;
            for (i = l + 1; i <= n - 1; i++)
            {
                d = a[i * n + l] / a[ll];
                for (j = l + 1; j <= n - 1; j++)
                {
                    u = i * n + j;
                    a[u] = a[u] - d * a[l * n + j];
                }
            }
        }
        return (k);
    }

#if __DRIVE_CODE__
      int main()
      { 
          double a[5][4]={ {1.0,2.0,3.0,4.0},
                                {5.0,6.0,7.0,8.0},
                                {9.0,10.0,11.0,12.0},
                                {13.0,14.0,15.0,16.0},
                                {17.0,18.0,19.0,20.0}};
          cout <<"RANK=" <<rank(&a[0][0],5,4) <<endl;
          return 0;
      }
#endif
}
